clear
matrix drop _all
capture log close

log using "suicide-1871-restat-final.log", replace

set matsize 800
set more off
set rmsg on


use "suicide-data-1871-restat.dta", clear

/****************************/
/* Define control variables */
/****************************/

global demog     "f_under15 f_over60 hhsize"
global urban     "f_urban"
global disab     "f_blind f_deaf f_dumb"
global econ      "perc_secBnC"
global educ      "f_rw"
global inter     "f_prot_f_urban"
global other     "f_fem f_jew f_ortsgeb f_pruss"  

#delimit;
global occ1882   "perc_occ1882_var_2 perc_occ1882_var_3 perc_occ1882_var_4 perc_occ1882_var_5 perc_occ1882_var_6 perc_occ1882_var_7 perc_occ1882_var_8 
                  perc_occ1882_var_9 perc_occ1882_var_10 perc_occ1882_var_11 perc_occ1882_var_12 perc_occ1882_var_13 perc_occ1882_var_14 perc_occ1882_var_15 
				  perc_occ1882_var_16 perc_occ1882_var_17 perc_occ1882_var_18_1 perc_occ1882_var_18_2 perc_occ1882_var_18_3 perc_occ1882_var_18_4 
				  perc_occ1882_var_19 perc_occ1882_var_20_1 perc_occ1882_var_20_2 perc_occ1882_var_20_3 perc_occ1882_var_21 perc_occ1882_var_22 
				  perc_occ1882_var_23_1 perc_occ1882_var_23_2 perc_occ1882_var_23_3 perc_occ1882_var_23_4 perc_occ1882_var_23_5 perc_occ1882_var_25";
#delimit cr				  


/*************************/
/* Figures 1 and 2: maps */
/*************************/

/***************************/
/* Figure A1: Scatter plot */
/***************************/

# delimit ;
graph twoway (lfit suic_1869_71_pc f_prot) (scatter suic_1869_71_pc f_prot, mcolor(black) msize(small)), 
xtitle(Share of Protestants) ytitle(Suicides (per 100,000 inhabitants)) graphregion(fcolor(white)) legend(off) 
ylabel(, format(%5.0fc)) xscale(range(-.01 1.01));
# delimit cr
graph export "scatter_suic_prot_1871.png", replace


/*************************************************/
/* Table 1: Descriptive Statistics, Prussia 1871 */
/*************************************************/

# delimit ;

global summvars "suic_1869_71_pc suic_1869_71_pd f_prot $demog $urban $econ $educ kmwittenberg f_fem f_jew f_ortsgeb f_pruss $disab facc_1869_71_pc facc_1869_71_pd lat_rad lon_rad zupreussen";

# delimit cr

sum $summvars


/******************************************************/
/* Table 2: Protestantism and Suicide in Prussia 1871 */
/******************************************************/

/* OLS */
/*******/
regress suic_1869_71_pc f_prot, robust
regress suic_1869_71_pc f_prot $demog, robust
regress suic_1869_71_pc f_prot $demog $urban, robust
regress suic_1869_71_pc f_prot $demog $urban $econ, robust
regress suic_1869_71_pc f_prot $demog $urban $econ $educ, robust
xi: regress suic_1869_71_pc f_prot $demog $urban $econ $educ i.regbez, robust
regress suic_1869_71_pd f_prot $demog $urban $econ $educ, robust


/*************************************************************************/
/* Table 3: Instrumental-Variable Estimates using Distance to Wittenberg */
/*************************************************************************/

/* 1st stage */
/*-----------*/
regress f_prot kmwittenberg, robust
regress f_prot kmwittenberg $demog $urban, robust
regress f_prot kmwittenberg $demog $urban $econ, robust
regress f_prot kmwittenberg $demog $urban $econ $educ, robust

/* 2nd stage */
/*-----------*/
ivreg2 suic_1869_71_pc (f_prot=kmwitt), robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ, robust


/***********************************************************************************************/
/* Table A2: Robustness to Conditioning on Potential Non-religious Outcomes of the Reformation */
/***********************************************************************************************/

ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ f_enrol1864_elem f_enrol1864_sec f_uni1871_enrol, cluster(kreiskey1864)
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ lnyteacher1886 tax1877_inctaxpc share_poor1901, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ ratio_yteacher1886_daylabor1901, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ n_protest_pre1871, robust


/**********************************************/
/* Table A3: Robustness to Additional Factors */
/**********************************************/

xi i.zupreussen

ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ $other, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ $other $disab, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ $other $disab f_married, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ $other $disab elevation, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ lat_rad lon_rad latlon, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ _Izupreusse_*, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ if county_dprot_l2==1, robust


/**********************************************************/
/* Table A4: Accounting for Unpleasant Weather Conditions */
/**********************************************************/

reg rain_n kmwitt, robust
reg temp_n kmwitt, robust
regress f_prot kmwitt $demog $urban $econ $educ rain_n temp_n, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ rain_n temp_n, robust


/**************************************/
/* Table A5: Alternative Death Causes */
/**************************************/

replace dead_1869_71_pc = dead_1869_71_pc/10

ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ facc_1869_71_pc, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ dead_1869_71_pc, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $educ $occ1882, robust
ivreg2 facc_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ eastprussia, robust
ivreg2 dead_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ, robust
ivreg2 suic_1869_71_pd (f_prot=kmwitt) $demog $urban $econ $educ, robust


/*************************************/
/* Table A6: Religious Concentration */
/*************************************/

ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ herfindahl, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ if f_prot<0.02 | f_prot>0.98, robust
ivreg2 suic_1869_71_pc (f_prot=kmwitt) $demog $urban $econ $educ if f_prot<0.005 | f_prot>0.995, robust


/******************************************************/
/* Table A7: descriptive statistics from Hilse (1871) */
/******************************************************/

/************************************************************/
/* Tables A8 and A9 and Table 4: see suicide-1816-restat.do */
/************************************************************/

/**************************************************/
/* Table 5: see suicide-1816-1871-panel-restat.do */
/**************************************************/

/**********************************************************************************/
/* Table 6: Discriminating between Sociological and Theological Explanations 1871 */
/**********************************************************************************/

global demog     "f_under15 f_over60 hhsize"
global urban     "f_urban"
global disab     "f_blind f_deaf f_dumb"
global econ      "perc_secBnC"
global educ      "f_rw"

global soci_theovars "f_comm f_prot q1_f_comm f_prot_q1_f_comm q4_f_urban f_prot_q4_f_urban"

for any f_urban: sum X, d \ gen q1_X = (X<r(p25)) \ sum X, d \ gen q2_X = (X>=r(p25) & X<r(p50)) \ sum X, d \ gen q3_X = (X>=r(p50) & X<r(p75)) \ sum X, d \ gen q4_X = (X>=r(p75))
for any  q1_f_urban q2_f_urban q3_f_urban q4_f_urban: gen f_prot_X = f_prot * X

for any f_comm: sum X, d \ gen q1_X = (X<r(p25)) if f_comm!=. \ sum X, d \ gen q2_X = (X>=r(p25) & X<r(p50)) if f_comm!=. \ sum X, d \ gen q3_X = (X>=r(p50) & X<r(p75)) if f_comm!=. \ sum X, d \ gen q4_X = (X>=r(p75)) if f_comm!=.
for any q1_f_comm q2_f_comm q3_f_comm q4_f_comm: gen f_prot_X = f_prot * X

regress suic_1869_71_pc f_comm $demog $urban $econ $educ if f_prot>.98, vce(cluster clus_f_comm)
regress suic_1869_71_pc f_prot q1_f_comm f_prot_q1_f_comm $demog $urban $econ $educ , vce(cluster clus_f_comm)
regress suic_1869_71_pc f_prot q4_f_urban f_prot_q4_f_urban $demog $econ $educ if f_comm!=., vce(cluster clus_f_comm)
regress suic_1869_71_pc f_prot q1_f_comm f_prot_q1_f_comm q4_f_urban f_prot_q4_f_urban $demog $urban $econ $educ , vce(cluster clus_f_comm)


log close

